arXiv:1505.01785v3 [nlin.PS] 1 May 2016 


Dam break problem for the focusing nonlinear Schrodinger equation and the 

generation of rogue waves 

G.A. El^, E.G. Khamis^’^, and A. Tovbis"^ 

^ Department of Mathematical Sciences, Loughborough University, Loughborough LEll 3TU, UK 
^ Instituto de Fisica, Universidade de Sao Paulo, 05508-090 Sao Paulo, Brazil 
^ Center for Weather Forecasting and Climate Studies-CPTFC, 

National Institute for Space Research (INPF), Cachoeira Paulista, Sao Paulo, Brazil 
^ Department of Mathematics, University of Central Florida, USA 
(Dated: May 3, 2016) 

We propose a novel, analytically tractable, scenario of the rogue wave formation in the framework 
of the small-dispersion focusing nonlinear Schrodinger (NLS) equation with the initial condition in 
the form of a rectangular barrier (a “box”). We use the Whitham modulation theory combined 
with the nonlinear steepest descent for the semi-classical inverse scattering transform, to describe 
the evolution and interaction of two counter-propagating nonlinear wave trains — the dispersive dam 
break flows — generated in the NLS box problem. We show that the interaction dynamics results 
in the emergence of modulated large-amplitude quasi-periodic breather lattices whose amplitude 
profiles are closely approximated by the Akhmediev and Peregrine breathers within certain space- 
time domain. Our semi-classical analytical results are shown to be in excellent agreement with the 
results of direct numerical simulations of the small-dispersion focusing NLS equation. 

PACS numbers: 


I. INTRODUCTION 

There has been much interest over the last two decades in the modelling rogue wave formation in the framework of 
the one-dimensional focusing nonlinear Schrodinger (NLS) equation, 

ieipt + = 0, (1) 

where if is the complex wave envelope and e is a free parameter defining the modulation scale; all variables are 
dimensionless. Rogue waves represent the waves of unusually large amplitude whose appearance statistics 

deviates from the Gaussian distribution. The conventional “rogue wave criterion” is > 2, where \ip\s is the 

significant wave height computed as the average wave height of the largest 1/3 of waves (see e.g. [1], [2], [3]). For the 
random wave held with Gaussian statistics one has \ip\g = 2|^/oPj where i/jq is the background (mean held) amplitude, 
leading to the criterion \if\'^ > 8|'0oP- 

As is clear from the above discussion, the description of rogue waves is a two-sided problem: it is concerned both 
with the dynamics of their formation and evolution, and with the statistics of their occurrence. In this paper we 
shall be concerned only with certain dynamical aspects of the rogue wave generation. The NLS equation (1) has a 
number of relatively simple exact solutions which are widely considered as prototypes of rogue waves, the principal 
representatives being the Akhmediev, the Kuznetsov-Ma and the Peregrine breathers (see, e.g., [4]). The main physical 
contexts for rogue waves are oceanography and nonlinear optics (see [1] , [5] , [6] and references therein). 

Finding controllable ways to excite rogue waves has been one of the central topics of the “deterministic” rogue wave 
research (see e.g. [9] and references therein). Various mechanisms have been proposed in the framework of the NLS 
equation and its generalisations (see e.g. [2], [7], [8] and references therein). Many of them relate the rogue wave 
appearance to the development of modulational instability of the plane wave due to small perturbations (see, e.g., 
[10], [3], [11]) or large-scale initial modulations [12]. Other proposed mechanisms involve nonlinear wave interactions: 
e.g., the interactions of individual solitons [13] or the interaction of solitons with the plane wave [14]. One should note, 
however, that, while there have been many papers developing the methods for finding particular rogue wave solutions 
(Darboux transformation, 9-method and oth. — see, e.g., [15], [14] and references therein), an analytical description 
of their formation from reasonably generic inital data remains a challenging, and to a large degree unsolved problem. 
In most cases numerical simulations remain the only available resort. 

One can distinguish two contrasting general classes of initial-boundary value problems associated with equation 
(1). The first one is concerned with the evolution of rapidly decaying potentials. The result of such an evolution 
generically represents a combination of fundamental solitons and some dispersive radiation with no rogue waves at the 
output. A comprehensive description of this process is achieved in the framework of the Inverse Scattering Transform 
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(1ST) [16]. The second class of problems deals with the NTS equation with non-decaying boundary conditions, and 
is much less explored analytically. One of the most interesting and physically relevant problems of this kind is the 
evolution of various perturbations of the plane wave (the “condensate”), 

^ ( 2 ) 

where the amplitude g > 0. Small harmonic perturbations ^ Qi(.kx-u]t) 1 - 2 ^ satisfy the dispersion relation uj{k) = 
±|fc[(efc)^ — 4 ( 7 ^]^/^, which implies modulational instability for sufficiently long waves with the wavenumbers k < 
kc = 2g/e. The description of the nonlinear stage of the development of modulation instability has been the subject 
of many papers including some very recent developments [11], [3], [17], [18]. It has been proposed in [3] that the 
long-time evolution of this process, in the case of random initial conditions, results in the establishment of a complex, 
globally incoherent strongly nonlinear wave state which can be associated with “integrable turbulence”, the notion 
introduced by Zakharov in [19]. It has also been shown that the rogue wave formation plays an important role in the 
characterisation of the early stage of the integrable turbulence development from a randomly perturbed plane wave 
[3], but also has a noticeable effect on the power spectrum of the established integrable turbulence developing from 
random initial conditions which are not small perturbations of a plane wave [ 20 ]. 

In this paper, we consider the problem that combines some of the key features of the both above contrasting funda¬ 
mental mathematical and physical set-ups (NTS with decaying vs. non-decaying boundary conditions). Specifically, 
we consider the evolution of a large-scale (compared to the medium’s coherence length) decaying data in the form of 
a rectangular barrier (a ‘box’ ) of finite height q > 0 and the length 2L: 


tp{x,0) 


q for jxj < L, 
0 for jxj > L. 


( 3 ) 


In the dimensionless variables of (1) we assume that L = 0(1), and the dispersion parameter £ <C 1. We shall refer 
to the problem (1), (3) as the small-dispersion NTS box problem. It is expected that the evolution (1), (3) at times 
t will model some features of the nonlinear stage of the development of modulational instability, in particular, 

the formation of rogue waves. 

The initial evolution of the box data for the small-dispersion focusing NLS equation can be viewed as a combina¬ 
tion of two “dam break” problems, which are known to lead to the instantaneous formation single-phase nonlinear 
modulated wave trains regularising the discontinuities at the opposite edges of the initial profile [21], [22]. These 
wave trains have the structure similar to dispersive shock waves (DSWs) [23], or undular bores, extensively studied 
in the defocusing (hyperbolic) NLS theory [24-27]. There are, however, important differences, which we emphasise 
by using the term dispersive dam break flow rather than DSW. The reason for using this term in the context of the 
focusing NLS equation is that the regularisation of discontinuous data via a single-phase modulated wave train in 
the focusing NLS occurs only if the upstream constant state is the vacuum, which the key feature of the classical 
shallow-water dam break problem [29] . However, the shallow-water “dry bottom” dam break problem does not involve 
the formation of a shock [29] so the regularising dispersive wave trains in the NLS box problem do not have classical 
shock counterparts in the dispersionless hyperbolic case. In contrast, the “genuine” focusing DSW analog is expected 
to have a multi-phase structure [30]. 

The dispersive dam break flows regularising the box data (3) expand inside the interval —L < x < L and, after 
a certain moment of time, to = ^/( 2 '\/ 2 'z), start to interact resulting in the formation of a region filled with a two- 
phase, x-t quasi-periodic wave, which we term the breather lattice due to the characteristic shape of the individual 
oscillations resembling standard breather solutions of the NLS equation. Indeed, we show that the wave form (the 
amplitude profile) of the oscillations in this lattice at each given t is quite well approximated by that of the Akhmediev 
breather with the spatial period depending on the value of t. Towards the end of the interaction region the spatial 
and temporal periods of the breather lattice increase so that locally, the oscillations are closely approximated by 
the Peregrine solitons with the amplitude 3q. Our numerical simulations show that further in time, the evolution 
leads to the generation of multiphase regions in the x-t plane with the number of oscillatory phases g (the “genus” 
of the solution) at any particular point —L < x < L growing with time, the expected asymptotic behaviour being 
g ~ t for t ^ 1 . Assuming the existence of the long-time asymptotics for the solution 'i/;(x,t,e) one can associate 
it with the ‘breather gas’ (see [31-34] for the description of the counterpart soliton gas in the KdV theory). One of 
the numerically observed features of the regions with (7 > 4 is the presence of the higher-order rogue waves with the 
maximum wave height Aq < < 5g. 

We now outline the analytical approach adopted in this paper. Although the NLS box problem admits exact 
analytical description via the 1ST [16], such a description becomes not feasible in practical terms when e 1 as the 
number of (global) degrees of freedom in the 1ST solution is then 0{e~^) 3> 1. In that case, the natural analytical 
framework is the semi-classical approximation, which enables one to asymptotically reduce the complicated exact 1ST 
construction to a more manageable description of a modulated multi-phase NLS solution approximating the exact 
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solution and involving just a few (local) degrees of freedom. Unlike in the exact solution, the 1ST spectrum of the 
approximate, semi-classical solution consists of finite number of bands which slowly evolve in space-time. There are 
two complementary mathematical approaches to the construction of such slowly modulated multiple-scale solutions 
to integrable nonlinear dispersive equations. The first one is based on the Whitham averaging procedure [28, 29] 
leading to a system of quasilinear equations governing the slow evolution of the endpoints of spectral bands [35-37]. 
For the focusing NTS, the Whitham modulation system is elliptic, implying modulational instability of the underlying 
nolinear periodic wave [36]. (One should be clear that, for the NTS equation, which is a modulation equation itself, the 
Whitham equations describe the ‘super-modulations’, i.e. the modulations of the nonlinear periodic or quasi-periodic 
solutions. See [38], [61] for the discussion of the relation between the NTS equation and the Whitham theory.) 

While the type (hyperbolic vs. elliptic) of the Whitham system yields the essential information about stabil¬ 
ity/instability of the underlying periodic wave [29], the modulation solution provides the detailed information about 
the evolution of a slowly modulated wave train. Vast majority of papers on the integration of the Whitham equations 
deal with the hyperbolic case, most notably in the DSW theory (see [23] and references therein). In contrast, the 
solutions of elliptic Whitham equations are far less explored, especially in the context of applications. The existing 
applied results are restricted to the simplest self-similar solutions in the single-phase case (see e.g. [39], [40], [41], [21], 

[42] ). 

The second approach to the semi-classical NTS equation is the nonlinear steepest descent method by Deift and Zhou 

[43] involving the Riemann-Hilbert problem (RHP) formulation of the 1ST. As a matter of fact, the two approaches 
are consistent, the modulation solution directly arising as part of the more general and mathematically rigorous (but 
also more technically involved) RHP analysis. In this paper, we use an appropriate combination of the two methods 
to construct a relatively simple analytic solution describing the x-t evolution of the approximate, semi-classical, 1ST 
specrum in the NTS box problem in the region of interaction of two dispersive dam break flows. This solution describes 
slow modulations of the two-phase breather lattice and enables us to predict the formation of rogue waves. We note 
that rigorous RHP analysis of the initial stage of the box evolution involving genus zero and genus one solutions was 
done in the recent work [22]. Part of the results obtained in [22] appear in the earlier papers [39], [21], where they 
were derived via the Whitham modulation theory. 

The semi-classical focusing NTS has only started to be explored as an analytical platform for the rogue wave 
research. We mention two recent papers [12] and [44] demonstrating the relevance of the semi-classical NTS scaling 
to the deep water ocean dynamics and the experimentally achievable configurations of Bose-Einstein condensates 
(BECs). Both above-mentioned works use the rigorous (RHP) mathematical results of [46] establishing the role of the 
Peregrine solitons in the regularisation of the focussing gradient catastrophe during the semi-classical NTS evolution 
of a certain family of analytic initial data that includes ijj{x,0) = sech(a;). Our paper continues this emerging line 
of research by putting forward and studying a novel scenario of the rogue wave formation via the interaction of two 
modulationally stable nonlinear wavetrains. Integrability of the NTS equation (1) and the semi-classical approximation 
play the key roles in the mathematical description of the proposed scenario. The proposed mechanism of the rogue 
wave formation can be realised in fibre optics experiments. Our analytical results are favourably compared with 
direct numerical simulations of the small-dispersion NTS box problem. In this regard we note that, although the 
semi-classical focusing NTS has been the subject of many numerical investigations (see e.g. [47-50] and more recent 
papers [51, 52] and references therein), we are aware only of a few examples, such as [53], [46] of the previous work 
undertaking quantitative comparison of the semi-classical analytical solutions for g > 1 with the direct numerical 
simulations of the NTS. 

The paper is organised as follows. In Section 2 we present an account of the necessary results from finite-gap theory 
of the focusing NTS equation and the associated Whitham modulation theory. Along with the well known results, 
this section contains a new general compact representation (21) for the characteristic speeds of the multiphased 
NTS-Whitham modulation equations. In Section 3 the modulation solution of the dam break problem is constructed 
following earlier analytical results of [39], [21], and then compared with numerical simulations. We show that, despite 
generic modulational instability in the system, this solution has the enhanced stability properties due to vanishing of 
the imaginary parts of the nonlinear characteristic speeds. Section 4 is central and is devoted to the analysis of the 
dispersive dam break flow interaction in the semi-classical NTS box problem. Using a combination of the Whitham 
modulation theory, the generalised hodograph transform [54] and elements of the RHP techniques we construct and 
analyse the modulation solution describing the interaction region, and then use it for the prediction of the rogue 
wave appearance. In Section 5 we numerically consider effects of small perturbations on the qualitative structure 
of the small-dispersion NLS box problem solution. We consider both perturbation of the initial conditions and 
perturnbations of the equation itself. In Section 6 we draw conclusions from our study and outline directions of future 
research. Appendix A contains an outline of the RHP calculations used in the construction of the modulation solution, 
and Appendix B is a brief description of the numerical method used for the solution of the NLS equation with small 
dispersion parameter. 
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II. MULTI-PHASE SOLUTIONS, ROGUE WAVES AND MODULATION EQUATIONS 


As is widely appreciated (see e.g. [55], [5]), one of the natural mathematical frameworks for the description of the 
development of modulational instability is the finite-gap theory [56] which is a counterpart of the 1ST for the NLS 
periodic problem [57]. The general finite-gap (or, better, finite-band) solution of (1) is given by 


i^g = q 


0g(a[/M/eU^ 

Qg{x/e,t/e;v'l) 


(4) 


where Qg is the Riemann theta-function associated with the hyperelliptic Riemann surface Tg of genus g given by 


Tlg{X-, a, a) = y^(A - q;o)(A - ao) • ■ ■ (A - ag)(A - ag), (5) 

where A is the complex spectral parameter. The branch points a. = («□) ■ ■ • j otg) and c.c. are the points of simple 
spectrum of the periodic non-self-adjoint Zakharov-Shabat scattering operator (see [56], [55], [5]). The phases 
in (4) are defined by the initial conditions. The plane wave solution (2) corresponds to the zero-genus spectral surface 
specified by Eq. (5) with = iq, i.e. TZo{X; a, a) = -y/(A — iq)(X + iq). Thus the spectral portrait of the plane wave 
is a vertical branch cut between the simple spectrum points ao = iq and do = —iq- (Note that, by associating a 
spectral portrait with a particular branchcut configuration we do not imply that the actual spectral bands (the spines 
in the terminology of [55], [5]) are necessarily located exactly along the branchcuts.) 

For g >1 the theta-solution (4) is a quasi-periodic function depending on g nontrivial oscillatory phases e~^r]j{x, t), 
so that V'g(- ■ • + 27r,.. •) = ■ • ■) for all j = 1 ,..., g. The phases are given by rjj = kjX -I- ujjt -|- 

j = 1, ... ,g. Here the (normalised by e) wavenumbers kj and the frequencies ojj are defined in terms of the branch 
points aj; and rfj are arbitrary initial phases. Also, associated with the ‘carrier’ plane wave is an extra, trivial phase 
770 = q^t. 

We present here the expressions for the wavenumbers kj and the frequencies ujj [36], [55], [5] which will be needed 
in what follows. 


kj = —‘iiri^Cj^i , 


ujj = —Airi 


g 

5 + ^k) Xj,l + Xj,2 , 

. k=0 


j = 


( 6 ) 


where Xj^k{oL^ a) are found from the system 



i^g-i 

'R.g{C,a,a) 


d-C — djk , 


j,k=l,...,g. 


(7) 


Here Stk is the Kronkecker symbol and % is a clockwise oriented loop around the branchcut connecting d/c and ak- 
Explicit expressions for Ki j for the case g = 2 can be found in Appendix A (see (A13)). 

For g = 1 the solution (4) is periodic and can be expressed in terms of the Jacobi elliptic functions. We introduce 
the notation aj = aj J- ibj. Then for the intensity we have (see e.g. [21]) 


IV’P = {bo + bi)^ - 4bobi sii^ (V(ao - ai)^ -h (&o + h)'^ {x-Ut + ^o)e , 

where the phase velocity U and the modulus 0 < m < 1 are given by 

(ao - Ao)(ai - di) 46o6i 


U = 4(ao -I- ao J- ai -I- ai) = oq J- oi, m = 


(ai - ao)(ao - ai) (oq - ai)^-|-(60 J-fei)^ 


( 8 ) 

(9) 


and ^0 is an arbitrary initial phase. For the wavenumber k of the cnoidal wave ( 8 ) we have 


Kirn) 


T^( x \/(ao - Qi)^ + (^0 + ^ 1 )^ , 
K [m) 


( 10 ) 


where K{m) is the complete elliptic integral of the first kind. Note that (10) can be obtained from the general 
representation ( 6 ), (7) by setting <7 = 1. The wave frequency w = kU then follows from the second expression ( 6 ), 
where = 0 . 
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In the harmonic limit, m = 0, the spectral portrait of the solution is a double point (either a\ = di or = do) 
on the real axis - this corresponds to a stable linear modulation of the plane wave ( 2 ); while the soliton limit, m = 1 , 
coresponds to two complex conjugate double points: ai = Oq and c.c. on the imaginary axis. This limit of (4) 
describes a fundamental soliton riding (or resting) on a zero background. 

The ‘classical’ prototypes of rogue waves (Akhmediev, Kuznetsov-Ma and Peregrine breathers) represent unsteady 
solitary wave modes on the finite background '0 = 9 and are described by special limits of the genus two (two-phase) 
solution (Eq. (4) with 5 = 2) [5]. The Akhmediev breather solution of (1) has the form (see e.g. [4]) 


cosh(fIte ^ — 2i0) — COS 0 cos(pa:£ 
cosh(nte“i) — cos(/)cos(pxs~^) 


( 11 ) 


where 


p = 2 sin 0 , and n = 2 sin( 20 ), ( 12 ) 

for 0 real. Solution (11) is periodic in space with the period P = 2tt£(p, and tends to the plane wave solution (2) in 
the limits t —>■ ±oo. The largest modulation occurs at t = 0 with the maximum envelope at x = 0, 

10^1 = 1-I-2 cos 0. (13) 

The Peregrine breather can be obtained from the Akhmediev breather by letting p —>■ 0. It is described by the 
rational solution of ( 1 ) 


0 P = 


4(1 -h 2it£-^) 

1 -|- + P) 


(14) 


which is localised both in space and time around x = 0,t = 0. It can also be obtained directly from the finite-band 
solution (4) with 5 = 2 by setting = 02 = ai = iq (and the c.c. expressions), i.e the spectral portrait of the 
Peregrine breather consists of two complex conjugate double points placed at the endpoints ±iq of the basic branchcut. 
The wave (14) has the maximum height |0p|max = Sg and represents a homoclinic solution starting from the plane 
wave ( 2 ) a.t t ^ —00 and returning to the same state at t —-l-oo (see [ 10 ] for the discussion of the special role of the 
Peregrine breather in the theory of rogue waves). 

The Madelung transformation 


0 = , 


(15) 


maps the NTS equation (1) to the dispersive hydrodynamics-like system with the negative classical pressure p = —p^/2. 


Ut + UUx - Px - 


pt + {pu)x = 0, 


( ^ 
V V 4p ) 


= 0 . 


(16) 


Here p > 0 and u are the density and velocity respectively, of the “fluid”. 

The prominent feature of the small-dispersion NLS evolution for decaying potentials is that, although it is globally 
described by the 1ST solution with a very large (~ ^ 1) number of degrees of freedom, the semi-classical 

asymptotics for t = 0(1) is locally (i.e. for Ax, At ^ e) described by the finite-band formula (4) with just a few 
degrees of freedom, but with slowly varying parameters. Specifically, the approximate solution (4) of the semi-classical 
NLS equation (1) is associated with the Riemann surface (5) of genus g = 0(1), whose slow (i.e. occurring on the 
scale much larger than e) spatiotemporal deformations are governed by the Whitham modulation equations for the 
branch points aj{x,t): 


{aj)t = vj‘^\a,a){aj)x, iaj)t = V^f\a, a){aj)x, j = 0,...,g. (17) 

Remarkably, aj are Riemann invariants of the Whitham modulation system (17), whose characteristic speeds 
Vj^\a,a) can be expressed in terms of Abelian differentials [36]. Another compact and physically insightful repre¬ 
sentations for 10-’s as nonlinear group velocities is (see [58], [59]) 

y{g) _ diOj ! dkj 
^ daj' dag ’ 


for any i = l,...,g. 


(18) 
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Formula (18) follows from the consideration of the system of g wave conservation equations 




(a, a) 


dx 


uJiia, OL 


j = 


( 19 ) 


as a consequence of the diagonal system (17). Equations (19) represent the consistency conditions in the formal 
averaging procedure [28], leading to the Whitham equations (17). In this procedure, the local conservation laws of 
the NLS equation (1) are averaged over the finite-band solutions (4) (see [36, 37]) respectively. On the other hand, 
equations (19) require that the fast phases s~^r]j{x,t) in the modulated finite-band solution (4) must be consistent 
with the generalised definitions of kj and iXj as the local wave number and local frequency respectively [29], [60], 


kj = {r]j)x, Wj = j = 1, • ■ • ,5- (20) 

Substitution of rjj = kjx + ujjt + in (20) yields the expressions for the initial phases ? 7 °, which are also subject to 
slow modulations and are expressed in terms of the spectrum branch points, ry°(a;,t) = Tj(q:,q:, ), j = 1,g [61]. 
As we shall show, the modulation phase functions Tj(a, a,) play the key role in the description of the rogue wave 
formation due to the interaction of the “regular” single-phase modulated wave trains. 

Using (18), (6), (7) one obtains the explicit compact representation for the characteristic speeds (from now on we 
omit the subscript g in TZg) 



= Re 



SLl ^fc,2 (C-4)-R(C) 

J2k=l (A-4W(C) 


( 21 ) 


where the parameters ^k ,2 are defined by (7). The derivation of the expression (21) for g = 2 is presented in 
Appendix A (see (A18)). 

The characteristic speeds (21) are generally complex-valued, i.e. the modulation system is elliptic implying modu- 
lational instability of the underlying nonlinear periodic wave [29] (see [62] for the historical account and [63] for recent 
advances in the mathematical understanding of the predictions of Whitham’s theory). However, it turns out that, 
for g > 1 some characteristic speeds (21) can undergo a degeneracy and assume real values (see (27), (26) below), so 
stable (or weakly unstable) configurations of nonlinear modulated waves in the focusing NLS dynamics are possible 
despite generic modulational instability. 

In the genus zero case, g = 0, the NLS modulation system (17) has the form 

(ao)i = (fao + 5 do)(ao)a;i (do)t = (fdo -f |Q!o)(Q;o)a: , (22) 

and is equivalent to the dispersionless limit of (1) 

pt + {pu)x = 0, ut + uux - px = 0, (23) 

where the Riemann invariants and characteristic speeds in (17) are expressed in terms of the hydrodynamic density 
and velocity as 


«o = -(2 +*V^). 


= fao + = -{u + iy/p). 


(24) 


One can see that the characteristics of (23) are complex unless p = 0 implying nonlinear modulational instability of 
the NLS equation (1) in the long-wave limit, which agrees with the linearised theory result for the plane wave (2). 
Since for p > 0 system (22) is elliptic, the initial-value problem for (22) is ill-posed for all but analytical initial data. 

For g = 1 the characterstic speeds (18) can be explicitly represented in terms of the complete elliptic integrals 
K(m) and E(m) of the first and the second kind respectively [37] 


= U + 


= u + 


(ao — ai){ao — ao) 


(ao - ai) + {ai - ao)E{m)/K{m) ’ 
{ai — ao)(ai — di) 

(ai - ao) -I- (ao - ai)E(m)/K{m)' 


(25) 


Here the modulus m and the phase velocity U are given by (9). Of particular interest are behaviours of the charac¬ 
teristic speeds (25) for m —>■ 0 (linear limit) and m —?> 1 (soliton limit). The linear limit can be achieved in one of the 
two ways (see (9)): either via ai —>■ di or via ao —>■ do. In the first case we have 


= 


— 2“° 2°^° ' 


u/i) = = 2ai 


52 

'^0 


a-i - ao 


ai —> ai : 


(26) 
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Similarly, 


7 2 o 1 

Qfo ^0 • ^ = ^0 ^ H- - —j ^1 = nOfi + -ai . (27) 

0-0 — 0-1 ^ z 

One can see that in both cases one complex conjugate pair of the characteristic velocities degenerates into a single 
real value while the other pair transforms into the pair of characteristic speeds of the genus zero system (22). 

In the soliton limit we have 


m —> 1 : 


ao = ai, =[/ = 2ai, 


(28) 


i.e. in this limit all the characteristic speeds are real, which is the modulation theory expression of stability of 
fundamental NLS solitons. 


III. DISPERSIVE DAM-BREAK FLOWS IN THE FOCUSING NLS EQUATION 

Within the modulation theory framework the initial conditions for the semi-classical NLS (1) are considered to 
be specified on the genus zero Riemann surface with added double points. The location of the double points in the 
complex A-plane is determined by the initial potential. The solution genus changes during the evolution implying 
the emergence of new oscillatory phases. In the spectral plane this process is the opening of the double points into 
spectral bands. The phase transition lines in the x-t plane, where the genus changes are often called the breaking 
curves. 



© 

Qo : 

ai=ai 

0 

f Qo : 

ai X 

Im K 


ReK ; 

ReK 


ai ^ 


Qo : 

; Qo : 


b) 




FIG. 1: Two contrasting generic types of the spectrum modification across the first breaking curve: a) transition (p = 0) —>■ 
(g = 2); b) transition {g = 0) ^ {g = 1) 


For a class of analytic initial data sufficiently rapidly decaying at infinity the semiclassical NLS evolution is initially 
defined on the genus zero Riemann surface and is governed by the elliptic system (22). This evolution leads to the 
onset of a gradient catastrophe, which is dispersively regularised via the generation of rapid, e-scaled, oscillations. 
These oscillations within certain neighbourhood of the gradient catastrophe point (xq, to) are asymptotically described 
by the finite-band solution (4) with g = 2 and slowly varying parameters aj [64] (see the proof of universality of such 
an asymptotic behaviour for analytic decaying potentials in [45], [46]). The spectrum modification across the first 
breaking curve in the above scenario consists in the opening of two Schwartz-symmetric bands (see Fig. la) so the 
solution in the region just above this curve describes a nonlinear two-phase wave which is prominently manifested 
in the appearance of the Peregrine soliton right beyond the point of the gradient catastrophe at a; = 0 [46]. The 
exact mechanism of the genus modification is described within the RHP approach in[64], [46] (see also [65]) but the 
necessity to introduce the genus two solution (rather than the genus one solution as it usually happens for “hyperbolic” 
dispersive systems like the defocusing NLS, and involves the DSW formation [25], [26]) can be qualitatively explained 
as follows. The genus zero system (22) is elliptic for all p > 0. One can readily see that this system does not support 
solutions involving variations of both variables p and u (or, equivalently og dg) along a single characteristic 
direction (in hyperbolic quasilinear systems such solutions are called simple waves). As a result, a generic gradient 
catastrophe in (22) necessarily occurs for a and a at the same point in the x-t plane, and so, its regularisation requires 
the introduction of two new, Schwartz-symmetric, spectral bands and thus, involves the change of the genus of the 
Riemann surface r(a:,t) from zero to two. 

There is another class of problems, where the phase transition across the first breaking curve occurs via the opening 
of a single spectral band emerging from the double point on the real axis (see Fig. lb), so the oscillatory solution 
{p{x, t), u(x, t)) asymptotically represents a modulated expanding single-phase wave train connecting the vacuum state 
(0,0) upstream with the plane wave (^^,0) downstream. This type of dynamics resembles the behavour of a DSW 
and is more in line with what one would expect in the hyperbolic case. As we shall see, this similarity is for a reason. 

















Consider the ‘dam-break’ problem for the focusing dispersive hydrodynamics (16): 


p(a:,0) 


> 0 for X < 0, 
0 for a: > 0, 


u{x, 0) = 0 . 


( 29 ) 


In the stable, defocusing, case the regularisation of the initial dam break (29) occurs via a smooth rarefaction wave 
and does not involve the generation of a nonlinear dispersive wave except for small (^ e) oscillations regularising the 
weak discontinuity at the rarefaction wave corner [25] ( see also [66]). In other words, the dam break problem solution 
to the defocusing NLS equation is asymptotically (e —0) equivalent to the classical dam break problem solution for 
shallow-water waves [29] . The dam break solution for the focusing NLS equation is of drastically different nature and 
involves dispersive regularisation via a finite-amplitude modulated single-phase wavetrain dehned inside an expanding 
transition region between two disparate states. In many respects this dynamic transition is analogous to a dispersive 
shock wave (DSW) [23] with the important difference that the ‘hyperbolic’ DSWs, similar to classical shocks, cannot 
provide transition from a vacuum state (see e.g. [29]). 

The modulation solution for the focusing dispersive dam break flow is constructed by noticing that one of the 
Riemann invariant c.c. pairs must be fixed to provide matching with the plane wave (2) downstream — this pair is 
ao = iq, do = The modulation solution for the second pair of the invariants ai{x,t), ai{x,t) must depend on 
s = x/t alone due to the scaling invariance of the problem (both initial conditions (29) and the modulation equations 
(17) are invariant with respect to the transformation x —> Cx, t —>■ Ct, where C is an arbitrary constant). As a result, 
the required modulation is given by a centred characteristic fan of the genus one Whitham equations (17) [39], [40], 
[ 21 ]: 


ao = iq, Re[Vi{iq,-iq,ai,ai)] 



lm[Vi{iq, -iq, ai, ai)] = 0. 


(30) 


(here and below in this section we omit the superscript for the characteristic speeds denoting the genus of the 
associated Riemann surface). Due to the symmetry in the expressions (25) for the characteristic speeds and the 
cnoidal solution (8) the modulation (30) could be obtained in terms of ao, ao using the characteristic speed Vo - in 
that case one would need to set constant the other pair of invariants: ai, di. 

We recall the notation aj = aj + ibj. Then, using (25) we obtain from (30) the explicit expressions 


af + (q-bif 
af-bf + g2 


E{m) 
K{m) ’ 


ao =0, bo = q, 

4g6i 

a? + (g + 6i)^’ 

X ^ q^ -bl 

--=2ai + ^ - 

t ai 


(31) 


Since the real part oi of the Riemann invariant enters (31) only as a square, the solution (31) defines two possible 
modulations of the cnoidal wave (8) corresponding to the right- and left-propagating dispersive dam break flows. 
To single out the modulation corresponding to the particular initial-value problem (29) it is instructive to represent 
solution (31) in terms of a single parameter 0 < m < 1 [39]: 


ai = ± 


2q 


mii{m 


■\/(l — m)[gL^(m) -I- to — 1], 


bi = 


mfj.(rn) 


[(2 — m)n{m) — 2(1 — m)], 


(32) 


- = ^ 

t m^{m) 


^/(i 


m){gt?‘{m) + m 


— / (2-TO)^(m)-2(l-m) \ 

\ lE{m)+m—l J ’ 


(33) 


where /i(m) = E{m)/K{m). Now one has to choose the lower sign in (32), (33) to provide correct matching with 
the plane wave (2) downstream, corresponding to the initial condition (29). The upper sign corresponds to the dam 
break flow of the opposite orientation, i.e. propagating upstream. 

It follows from (32) that for the chosen (leftward) propagation direction, 


m —^ 0 bi =0, ai = —g/-\/2; 


and 


1 oi = 0, bi = q. 


(34) 


Then, considering solution (31) in the limits m —> 0 and m —>■ 1 we obtain the speeds of the trailing (leftmost) and 
leading (rightmost) edges to be —2\/2qt and 0 respectively. Thus, the modulated wave train defined by (8), (31) 
is confined to the expanding region —2y/2qt < x < 0, where the modulus gradually varies from m = 0 (harmonic. 
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leftmost, edge) to m = 1 (soliton, rightmost, edge). The wave amplitude A = 2y/biq varies from A = 0 at the 
harmonic edge to A = 2q at the soliton edge where the wave form is described by the formula (assuming zero phase 
shift ^0 in ( 8 )) 


ijjs = 2,q sech{2qx/ . (35) 

At the harmonic edge x = —2\/2qt the wavenumber (10) assumes the value fco = k{m = 0) = \f%qle > 2qle (stable 
modulation). We note that the harmonic edge velocity s_ coincides with the linear group velocity of the left-going 
wave ui'ik) = —((efc)^ — 2 g^)((efc)^ — evaluated at fc = k^. 

Remarkably, solution (31) is essentially a simple-wave solution of the modulation system (17) with g — 1, for which 
two of the Riemann invariants (aojdo) are constant, and the remaining two (ai,di) vary along the same double 
characteristics family V 2 = x/t. As we mentioned earlier, the genus zero system ( 22 ) does not support solutions of 
this type. 

Comparison of the modulated wave train described by ( 8 ), (31) with the numerical simulation of the dam break 
problem for the focusing NTS equation is presented in Fig. 2 and shows excellent agreement. We note that an 
adjustment in the position of the wavetrain was required due to the use of smoothed initial step data in the numerical 
simulations. Such a long-time persistence of the self-similar modulation dynamics in the smoothed step problem is 
well established for the stable, defocusing, case (see [67]), but is not obvious and quite remarkable in the present, 
focusing, case, and deserves further analytical justification. 



X 


FIG. 2: (Color online) Dispersive regularisation of the dam break (29) in the focusing NLS eq. (16): numerical (solid line) vs 
analytical (modulation theory, dashed line) solution for |'i/)| = ^/p. The parameter values used are: q — 1, e = 1/33. 

We now discuss stability of the obtained solution of the focusing dispersive dam break problem. An important 
observation is that the constructed solution ( 8 ), (31) essentially describes a modulated soliton train. Indeed, the 
dependence of the modulus m{x/t) in this solution shown in Fig. 3a exhibits the values of m close to unity almost 
over the entire oscillations zone, which means that the dispersive dam break flow in the focusing NLS is dominated by 
solitons with the amplitude close to 2q. We note that the approximation of a dispersive dam break flow in a focusing 
medium by a modulated train of solitary waves was successfully used in [ 68 ]. The transition from m « 1 in the bulk 
of the wave train to to = 0 at the harmonic edge occurs within a narrow (in x/t-coordinate) dynamic region, where 
the new oscillations are generated and then quickly transform into solitons. 

Despite the partial saturation of the modulational instability in the dispersive dam break flow due to the vanishing 
of the imaginary part of the characteristic speed Vi for the solution (31), the wave train described by this solution is 
still subject to the instability implied by the nonzero imaginary part for the second pair of the characteristic speeds 
Vq,Vq associated with the constant Riemann invariants Oq = * 9 ) do = —iq. This instability, however, has a weak 
effect on the dispersive dam break flow as such due to the just mentioned fact that the major part of the wave train 
is dominated by solitons, which are modulationally stable. To quantify the effect of dispersive saturation we compute 
the imaginary part of (we restore the upper index here to distinguish this speed from the value (24) in the 
genus zero region). Using (25) and (31) we obtain 

_ ‘^qg{m){q^ - bf) _ 

[ai(l - /x(to ))]2 -b [(g - bi) + (g -b 6 i)/x(to )]2 ’ 


7 = Im[V[|^^] 


(36) 
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where, we recall, ^{m) = E{m)/K{m) (see (31)) and ai(m), bi{m) and m{x,t) are given by (32), (33). The value of 
7 can be viewed as a growth rate of nonlinear mode associated with the spectral pair ag, do- 

The dependence (36) of 7 on t for a fixed ccq < 0 is shown in Fig. 3b. One can see at each point —2^/2qt < xg < 0 
the value of 7 decreases to zero with time, thus confirming our conclusion about stability of the constructed solution. 




FIG. 3: Stability of the dispersive dam break flow for the focusing NLS equation, a) Behaviour of the elliptic modulus m as 
function of x/t in the modulation solution (33) implying the dominance of solitons in the bulk of the wave train; b) The decay 
of the growth rate 7 = of the unstable nonlinear mode associated with the pair ag = *?, do = —iq of the spectrum branch 

points. The plot shows the function 7 ( 1 ) at a fixed point x = xg = —0.1 inside the wave train described by the modulation 
solution (31). 


While the initial dam break undergoes rapid dispersive saturation, the upstream uniform plane wave state (2) for 
X < —2\/2qt remains modulationally unstable with respect to long-wave (k < 2 g/e) perturbations, and such small 
perturbations (a noise) are inevitably present in any physical system or numerical simulations. As a result, this 
instability imposes restrictions on the ‘natural’ lifespan of the described single-phase coherent structure of the ‘ideal’ 
focusing dispersive dam break flow. Let the typical amplitude of the noise be 5 <C 1. Then the linear theory prediction 
for the characteristic time of the development of the fastest growing mode with k = V^qfe is 


t 


m 


e 1 1 
2q‘‘ 0 


(37) 


which gives an estimate for the lifetime of the focusing dispersive dam break flow. 

In conclusion of this section we note that the focusing dispersive dam break problem was studied experimentally in 
[69] in the context of diffraction from an edge in a self-focusing medium. The authors of [69] observed the expanding 
nonlinear oscillatory regularisation of a discontinuous intensity profile, qualitatively similar to that described by 
solution ( 8 ), (31). To suppress the modulational instability of the background state in [69] nonlocality was used as 
suggested by previous theoretical studies [70], [ 68 ]. Solution ( 8 ), (31) was also used in [42] for the modelling of the 
matter-wave bright soliton generation at the sharp edge of density distribution in a BEC. 


IV. INTERACTION OF FOCUSING DISPERSIVE DAM BREAK FLOWS AND THE GENERATION 

OF ROGUE WAVES 

As was pointed out in Section II, the prototypical rogue wave solutions (soitons on finite background) are associated 
with the degenerate genus two NLS dynamics. This suggests that one can expect the rogue wave formation in the 
processes involving the interaction of two “regular”, single-phase waves {g = 1). Indeed, the “elementary” rogue wave 
events during individual soliton collisions were observed in numerical simulations [13] (see also [ 8 ]). Here we consider 
a more general scenario where rogue waves (not necessarily exact Akhmediev or Peregrine breathers) are formed in 
the interaction of two single-phase dispersive wave trains. 

We consider the NLS equation (16) with £ ^ 1 and initial conditions in the form of a rectangular barrier for the 
intensity with zero initial velocity. 


p{x,0) 


q^ for \x\ < L, 
0 for ja;j > L; 


u(x, 0 ) = 0 . 


(38) 
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We shall refer to (16), (38) as the NLS box problem. The rigorous semi-classical asymptotics of the NLS box problem 
solution were calculated in the recent paper [ 22 ] for the initial stage of the evolution involving solutions with 5 = 0 , 1 . 
The analysis in [22] was performed using the steepest descent for the oscillatory Riemann-Hilbert problem (RHP) [43] 
associated with the semi-classical NLS equation ([71], [64]) and, in particular, yielded the modulation solution ( 8 ), 
(31). In what follows we shall use an appropriate combination of the Whitham theory and the RHP techniques to 
construct a relatively simple exact modulation solution describing the dispersive dam break flow interaction [g = 2). 
This solution will then be used to predict the rogue wave formation. 


A. Before interaction {g — 0,1) 

The semi-classical evolution (1), (38) starts at t = 0 with the instantaneous formation of two dispersive dam break 
flows at the discontinuity points x = ±L of the initial profile. The corresponding modulation solution consists of two 
centred fans emanating from x = ±L and described by the formulae (see (30)): 

X + Re[P/^^(f(7, —iq, ai^ai)\t = ±L , 

(i) - 

The upper sign solution (39) is defined for 0 < a: < L and the lower sign for —L < a; < 0. Also one uses ai = 3fiai < 0 
in (39) for 0 < a: < L and oi > 0 for —L < a; < 0 - see (32). We intentionally represented solution (39) in the 
‘hodograph’ form to elucidate the natural matching of (39) with the modulation solution in the interaction region, 
which is of our primary interest. Explicit expressions for the modulation solution (39) in terms of elliptic integrals 
are obtained from (31) by replacing x with a: ± L. 

Earlier we introduced the notion of a breaking curve as the line in x-t plane separating regions described by solutions 
with different genus g. On the first breaking curve Ti separating the regions with g = 0 and g = 1 (see Fig. 6a) one 
has (see (34)): 


Ti : m = 0 , q;i==F<?/V 2 , 

where the minus sign applies to 0 < a; < L and the plus sign to —L < x < 0. Now, evaluating the limit m 
solution (39), we obtain the equation of the first breaking curve 


Ti : 


t = 


L- \x\ 

2V2q 


(40) 
0 of the 

(41) 


consisting of two symmetric parts. 



FIG. 4: (Color online) Collision of two counter-propagating focusing dispersive dam break flows results in the formation of 
rogue waves. Numerical simulation of the NLS box problem (16), (38) with e = 1/33, q = 1, L — 25/33. The rogue wave is 
formed at about t = 0.376. a) Density plot for the amplitude j'i/ij = y/o; b) Amplitude prohles for different t. 
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At X = Ot = to = Ll2\f2q the counter-propagating dispersive dam break flows collide, so the developed single-phase 
description (39) becomes invalid for t > to. In Fig. 4 the results of the numerical simulations of the box problem 
with L — 25/33, q = 1 and e = 1/33 are presented for t < to (the smoothed initial data were used in the numerical 
simulations to avoid unwanted numerical effects). The numerical method used in the simulations is briefly described in 
Appendix B, where also the effects of smoothing the initial data are discussed. In Fig. 4a the density plot is presented 
for \ij}{x,t)\ while Fig. 4b shows spatial profiles of |'0(3^)l at different moments of time. One can see from the bottom 
plot that the dispersive dam break flow collision leads to the rapid formation of a rogue wave with the maximum at 
a: = 0, the amplitude slightly less than 3 and the wave form typical of a breather: a tall central peak rapidly decaying 
to zero and two smaller “wings” at both sides. 

Remarkably, at the point of collision (x = 0,t = to) the amplitudes of both single-phase modulated wavetrains 
are very small and the modulation is stable (the modulation solution yields the zero amplitude and the wavenumber 
ko = \f%qle > fcc at this point). However, for t > to the interaction between these two stable, small-amplitude tails 
of the counter-propagating dispersive wavetrains gives rise to the development of modulational instability resulting in 
the rapid formation of rogue waves. 

In the next section, we shall study the development of this process using the modulation theory for g = 2, some 
results of the RHP analysis of the semi-classical NTS equation and direct numerical simulations. 


B. Interaction (g > 2) 

We first present the results of the numerical simulations of the NTS box problem with q = 1, L = 25/33 and 
two different values of e: 1/33 and 1/60. The respective x-t density plots for the amplitude \ip\ = are presented 
in Fig. 5a and Fig. 5b. One can see the regions with distinctly different behaviour of the solution in both plots. 
Remarkably, both simulations produce very similar macroscopic patterns differing apparently only in the number of 
the oscillations forming the pattern. This striking robustness of the macroscopic features of the dispersive dam break 
flow interaction despite the presence of modulational instability in the system, is a strong indication of the applicability 
of the limiting (e —> 0), semi-classical description and of the Whitham modulation equations in particular. Indeed, 
as we mentioned, the existence of the semi-classical limit was rigorously established in [22] for the initial stage of 
the evolution involving solutions with g = 0,1 but now we have some confidence in assuming the validity of the 
semi-classical description for the regions with g > 1. 



FIG. 5: (Color online) Density plot for in the focnsing NLS box problem with q = 1, L = 25/33. a) e = 1/33; b) 

£ = 1/60. 

We have used the image filtering software (an edge-detect filter is applied to detect contours and facilitate the 
sampling of points) to extract the boundaries between the regions with qualitatively different behaviours of the 
oscillations in the plots of Fig. 5. The normalised to q = 1, L = 1 diagram is shown Fig. 6a where Bezier curves of 
a few points extracted from the filtered image are used to produce the boundary lines. A minor adjustment of the 
image to the analytically available points (such as the collision point (0, l/2-\/2)) was necessary to compensate for the 
slight time delay present in the direct simulations in Fig. 5 due to the smoothing of the box edges. 

We now formulate the hypothesis about the structure of the x-t plane in the box problem shown in Fig. 6a . As 
our numerical simulations suggest, the interaction of two single-phase dispersive dam break flows is described by the 
modulated two-phase {g = 2) solution confined to a curved rhombus-like region. The genus change pattern from 
g = 1 to g = 2 across the first breaking curve is also consistent with the spectrum modification mechanism shown in 
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Fig. lb as it involves the emergence of one band from a double point on the real axis of the A-plane. The opposite 
signs of the real parts of the evolving Riemann invariants in the two dam break solutions (39) {g = 1) before the 
interaction, suggest that the spectral portait in the interaction region will have a qualitative form shown in Fig. 6 b, 
i.e. will consist of two slowly evolving bands 71 and 72 located at opposite sides of the fixed central band 70 (see also 
Appendix A). Indeed, in Section 4.2.1 below, we shall construct the modulation solution describing the slow evolution 
of this spectral portrait and show that it is consistent (i.e. provides the necessary matching) with the known “outer” 
solution (39) in the 5 = 1 region. 

The subsequent evolution leads to the formation of the regions of higher genera: 5 = 3, 5 = 4, etc. (see Fig. 6 b)). 
All the generated oscillations are confined to the original box interval —L < x < L. Remarkably, for any 0 < |a;| < L, 
each crossing of the breaking curve results in the genus increment by one, while the genus increases by two when 
crossing the breaking curves at a; = 0 . We note that the known scenario of the development of the semi-classical 
NLS solution for analytic, sufficiently rapidly decaying initial data (e.g. 'ijj{x,Q) = sech(a;)) involves only the genus 
increments by two across breaking curves [64], [71] (see Fig. 7). The striking difference between the two scenarios 
is the reflection of the two fundamentally different spectral mechanisms of the genus transformation involved. These 
are shown in Fig. 1 for the Hrst breaking curve but the principle remains for higher breaks as well. 




FIG. 6: a) (Color online) The structure of the x-t plane of the NLS box problem in the semiclassical limit suggested by the 
x-t density plots in Fig. 5. The diagram is normalised for g = 1, L = 1; b) The spectral portrait of the solution in the g = 2 
region;. 



FIG. 7: NLS equation (1) with e = 1/33 and the analytic initial condition ip = sechx. a) (Color online) Density plot for 
\ip{x,t)\; b) The structure of the x-t plane in the semi-classical limit suggested by a). 


Below we concentrate on the analysis of the dispersive dam break flow interaction occurring in the region with 
5 = 2. 
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1. g = 2: hodograph solution 


The generic spectral portrait of the finite-band NLS solution in the g = 2 region of the box problem is shown in 
Fig. 6a. The modulation is described by three complex conjugate pairs of Riemann invariants aj,aj, j = 0,1,2. 
Two of the Riemann invariants remain constant: ao = *9: do = —iq-, so we are left with just two varying pairs: 

j = 1,2. We note that, within the framework of the modulation theory, the constancy of aQ,ao 
follows from the matching, across the second breaking curve, of the genus two modulation solution with the genus 
one solutiion (39) for which these invariants are constant. At the same time, in the rigorous RHP analysis of [22], 
the branchpoints ao, do for the box potential are special, as they coincide with the endpoints of the spectral interval 
on the imaginary axis (this spectral interval is the locus of accumulation of the points of discrete spectrum, whose 
number is growing like 0{e~^)). The spectral interval, of course, is invariant under the NLS-evolution. Further, due 
to the symmetry x —)■ —x of the box problem the evolving bands must always be located at the oposite sides of the 
imaginary axis Re[A] = 0, hence at a; = 0 one has ai = —a 2 - It is worth mentioning that the opposite signs for Re[ai] 
and Re[a 2 ] most readily follow from the modulation theory as they are requuired by the matching, along the second 
breaking curve, of the modulation solution in the genus two region, with two genus one modulation solutions for the 
dispersive dam break flows propagating in opposite directions (towards each other). By construction, the varying 
Riemann invariants for these flows have real parts of opposite signs - see Section 4.1. The same result follows from 
the rigorous analysis of [22]. The solution for the moving invariants aj,dj, j = 1,2 can be found via the Tsarev 
generalized hodograph transform [54]. This method was originally developed for hyperbolic systems of hydrodynamic 
type but is equally applicable to elliptic systems. The Tsarev result in the application to our present problem can be 
formulated as follows. Any local non-constant solution of the modulation system (17) in the genus two region is given 
in an implicit form by the system of two algebraic equations with complex coefficients 


X + Vj{a,a)t = Wj{cy.,cy.), x+ V = Wj{a,a), j = l,2, 


(42) 


(o') 

where a = {iq,ai,a 2 ); the characteristic speeds Vj{a.,Q.) = Vj ’{a.,d.) are given by (18), (6) (or, equivalently, by 
(21)). The four unknown functions Wj, Wj, j = 1,2 satisfy the system of four linear partial differential equations 


dajWk 
Wk - Wj 


dc,,Vk 

Vk - V, 


and C.C.; 


j.k = 1 , 2 , 




(43) 


where = gfj. 

For (42), (43) to describe the interaction of two dispersive dam break flows in the box problem, equations (43) 
must be supplied with appropriate boundary conditions. These conditions follow from the requirement of continuous 
matching, on the second breaking curve t = T 2 {x), of the hodograph solution (42) with the known solution (39) 
(in which ai should be replaced by 02 for 0 < a: < L in the “plus” branch of the solution). Similar to the first 
breaking curve Ti, the curve T 2 is a free boundary, on which two of the Riemann invariants merge (cf. Fig. lb 
showing the prototypical spectrum modification across a breaking curve in the box problem). It can be shown that 
the corresponding merged velocities are real (see [65]) so the boundary T 2 represents a symmetric curve consisting 
of two real characteristics that carry over the constant values of the merged Riemann invariants brought by the two 
branches of Ti (see (41)) 


T 2 : ai = ai = —qj\j2 for — L < x < 0 and 02 = 02 = q/V2 for 0 < x < L . 


(44) 


The matching regularisation procedure for the the focusing NLS is in many respects analogous to that in the defocusing 
NLS theory (see [72]). We won’t describe it in much detail but just mention that, having found the solution Wj, Wj 
of the Tsarev equations, one then needs to verify that the resulting hodograph equations (42) are invertible, i.e. that 
they specify functions a{x,t), d{x,t) in a certain region of x-t plane. See [73] for the ‘hyperbolic’ counterpart of this 
construction arising in the description of the DSW interaction in the KdV modulation theory. 

We by-pass the outlined above direct matching regularisation procedure by taking advantage of the available 
mathematical results from the RHP analysis [77] of the semi-classical focusing NLS, and applying them to the genus 
two region in the box problem. More specifically, we shall express Tsarev’s Wj’s in terms of the “modulation phase 
shift” functions Tj(Q:,Q;) introduced in Section 11. These phase functions depend in a simple way on the the semi- 
classical scattering data for the box potential [22]. We shall then verify that the functions Wj generated by the phases 
Tj-. i) satisfy equations (43); and ii) provide, via the hodograph formulae (42), the required matching for aj(x,t) on 
the breaking curve. Below we present the formal derivation of the modulation solution along these lines. The outline 
of the self-contained rigorous RHP construction underlying this derivation and explaining the connection between the 
modulation solution and the semi-classical RHP can be found in Appendix A. 
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We start with recalling that, for the solution t; e) of the semi-classical NLS to have an asymptotic representation 
in the form of the modulated finite-band potential (4), (17) depending on g phases of the form e~^rjj = e~^(kjX + 
ujjt + rjj) (i.e. linear functions of x and t) one must require that the “initial phases” rj^ are functions of a., a so 
that the general kinematic conditions (20) are satisfied. To this end we introduce rjj = —Tj{a.,a.) and use the first 
condition (20), 


to obtain 


d{kjX + LUjt 

dx 



J = l,2, 


(45) 


dkj 

do-jYi 


X + 


dujj ^ dTj 


and C.C., j, m = 1, 2, 


(46) 


provided d^ctj ^ 0, ^ 0, j = 1,2. Here kj{a,a) and ujj{cy.,a) are defined by (6), and T^’s are yet to be found. 

The second condition (20) leads to the same set of equations (46). As we shall see, only half of the equations (46) are 
independent, so it is sufficient to consider either j = 1 or j = 2. We also note that equations (46) admit a compact 
and elegant representation in the form of the stationary phase conditions: 


9rij 

dotyri 


= 0 , 


dr], 

dotjn 


= 0, j,m = 1,2. 


(47) 


For given functions Tj(q:,q:) equations (46) fully define the modulations a{x,t),a.{x,t) (assuming invertibility of 
(46), which is not guaranteed a-priori). Comparing equations (46) with the hodograph solution (42) and using the 
representation (18) for the characteristic speeds Vj{a,a) in (42) one readily makes the identification 

da T 

Wm = a"' J and C.C., j,m = 1,2. (48) 

^^06 jjl 

Now we observe that formula (48) must yield the same function Wmi’Ct, ct) for both values of j. This is a consequence 
of the consistency of the genus two Whitham modulation system with two “extra” (wave number) conservation laws 
(19) (the same argument was used to establish the “nonlinear group velocity” representation (18) for the characteristic 
speeds of the Whitham modulation system). Thus, it is sufficient to consider only half of the equations (46). 

Within the RHP approach the phases are determined in terms of g-l-1 functions /^(C), k = 11,1,... ,g containing 
the information about the scattering data for the initial potential (see Appendix A and [65]). For g = 2 we have 


where 


- _ / fkiQpjiQdC 

TZiC) 




(49) 


Pj{X) = XTjpA -I- 3<j,2, (50) 

77,(A) = 77.2(A, a, a) is given by (5), and the coefficients Xj ,2 are determined by conditions (7) with g = 2. 
We have verified that functions wj defined by (42), (48), (49) satisfy Tsarev’s equations (43) for arbitrary fk{C), 
thus defining the general local solution to the genus two Whitham equations. We mention in passing that the 
quantities dHj = pj(X)dX/TZ{X, a, a) represent the normalised holomorphic differentials playing important role in 
the construction of finite-band solutions of the NLS equation (see e.g. [55]). These differentials also serve as the 
generating functions for solutions of the Whitham equations associated with KdV hierarchy [58], [59]. One can see 
from (49), (48) that they play essentially the same role in the NLS modulation theory. 

For the box potential (38) we have ao = *9 so 77.g(A) = R{X)i'{X), where 

i?(A) = yJiX- ai){X - di)... (A - Q;g)(A - dg), i^(A) = y^X^~+~^. (51) 

The functions /fc(A) in (49) can be inferred from the results of [22]: 

/o(A) = /i(A) = -2LX, MX) = -2LX + 4Lu{X). (52) 

Using (48), (42) one can verify that the necessary matching conditions with the genus one solution (39) on the second 
breaking curve are satisfied. In what follows we shall only need to explicitly verify these conditions at the point of 
the dispersive dam break flow collision x = 0,t = L/2\/2q, which is the common point for T, and T 2 . At this point 
one has a, = —qj^M and 02 = qlM2 (see (44)). 
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2. Rogue wave formation 


We shall call the modulated quasi-periodic wave in the dam break flow interaction [g = 2) region the modulated 
breather lattice, by analogy with the term ‘modulated soliton lattice’ used for the slowly varying finite-band solutions 
of the KdV equation [60]. We first present numerical solution for the amplitude \if’{x,t) \ in the interaction region. The 
spatial profiles shown in Fig. 8 represent the snapshots of the solution in Fig. 5a taken at the times corresponding to 
the temporal maxima of the breather lattice. 



FIG. 8: (Color online) Emergence and development of the modulated breather lattice {g — 2) due to the interaction of dispersive 
dam break flows {g = 1). Numerical solution \ip{x,t)\ of the NLS equation (1) with e = 1/33 and the box initial conditions 
(38) with q — 1, L = 25/33. See Fig. 5a for the corresponding density plot. 


We shall use the modulation solution obtained in the previous subsection to analyse the dam break flow interaction 
dynamics. For that, we shall look at the temporal behaviour of the modulations at a: = 0, which, due to the symmetry 
X —> —X of the initial data, allows for a significant simplification of the analytic expressions. As we shall see, the 
modulation solution at a; = 0 provides a major insight into the interaction dynamics in the entire genus two region. 
From (46) we obtain: 


a: = 0 : 


dujj ^ dTj 


and c.c., m,j = 1, 2. 


(53) 


As we have already mentioned, in view of the symmetry x —> —x of the box problem the solution at a: = 0 
must also exhibit the spectral “portrait” symmetric with respect to the imaginary axis Re[A] = 0 so that we have 
0 ^ 2 ( 0 ,t) = —di(0,t), and the coefficients entering the expressions for u;i _2 and Ti ^2 (see (6), (7) and (49), (50)) 
can be evaluated as (see Appendix A) 


a; = 0 : >^1 1 = —>f 2 1 = — uv 

4i 

where 

Q{z) = \/[{z - 6)2 -b a2][(z -I- 6)2 -|- a^], g{z) = \/- z"^ . (55) 

As follows from our previous analysis, only half of the equations in the system (53) are independent so it is sufficient 
to consider only j = 1 or j = 2 which leaves us with two c.c. equations. The symmetry 02 = at a; = 0 reduces the 
number of independent equations to just one. As a result, the hodograph modulation solution (53) can be represented 
in the form (see Appendix A for the outline of the calculation) 


zdz 


Q{z)g.{z) 


Xl,2 — X2,2 — ~7r- 
2^ 


dz 


Q{z)g{z) 


(54) 


L 

271 


dz 


Q{z)n{z 


{z — b) — ia dz 
\z + iaY Q{z) 


{z — b) — ia dz j t L f dz 

|z-|-iap Q{z)fj,{z) 1 2 27r y Q{z) 


-g 


-g 


(56) 
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where 


a{t) = a{t) + ib{t) = a2{0,t). (57) 

It is not difficult to verify that at the collision point t = L/2^/2q one has a = qjV^, b — 0 as required by the matching 
conditions (41) and (44). Thus, the obtained solution (56) indeed satisfies the matching conditions at the breaking 
curve Ti. 




FIG. 9: Evolution of the spectral branch point 02 = —hi = a + ib at x = 0 in the genus two interaction region - formula (56) 
with g = 1, L = 25/33, as in Figs. 5, 8. The evolution starts on the real axis, b = 0, a = l/\/2 at t = to = L/2\/2 « 0.27. One 
can see that b rapidly grows with time while a decreases, so both branch points Q 2 and ai approach the imaginary axis close 
to 6 = 1 at the end of the genus two region, a) Dependencies a{t) and b{t); (b) Trajectory of 0:2 = —hi in the complex plane. 


Separating real and imaginary parts in (56) we obtain a system of two equations for a{t) and b{t). We then solve 
the obtained system numerically, using the Broyden method [74]. This method is a generalization of the secant 
method, also known as a quasi-Newton method, to systems of nonlinear equations, and involves the Jacobian matrix 
of the system. Instead of the analytical evaluation of derivatives in (56), we constructed an approximation to the 
Jacobian matrix, which was updated at each iteration. The initial Jacobian can be set as the identity matrix or 
a finite difference approximation, like in the first-order forward difference method. The iterations stop when some 
tolerance value, e.g. 10“®, is reached. The resulting plots a{t), b{t) and b{a) are presented in Fig. 9. One can see that 
the imaginary part of a rapidly grows, while the real one decreases. 




FIG. 10: (Golor online) a) Dependence of the spatial period (the wavelength) on time in the genus two interaction region at 
a; = 0. Solid line - formula (58); dots - the average distance between the peaks of the breather lattice adjacent to the central 
peak at a; = 0 in the numerical solution of the NLS box problem (g = 1, L = 25/33, e = 1/33), at different moments of time; 
b) Comparison of the numerical solution for the amplitude \ip\ in the box problem at f = 0.676 (solid line) with the (shifted) 
amplitude profile of Akhmediev breather (11) with the spatial period 27r£/p = P(0.676) = 0.154 found from (58) (dashed line). 

As follows from (6), (54), the ^-normalised wavenumbers and frequencies are calculated at a: = 0 as fci = —^2 = 
—darlxTi 1 and wi = a ;2 = —47ri>!:i^2 respectively. Thus, in the vicinity of a; = 0 the wave can be viewed as periodic 
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FIG. 11: (Color online) a) Density plot for in the NLS box problem with q = 1, L = 25/33, e = 1/33 - the zoom in on 
the region of the rogue wave formation at f = 1.166; b) The rogue wave profile: solid line - the numerical solution of the box 
problem zoomed in near x = 0, t — 1.166; dashed line - the Peregrine breather (14) for q = 1. 


with the spatial P and temporal (‘breathing’) T periods slowly depending on time 


Pit) 


OO 

27re f zdz 

ki ^ J Q(z)/r(z)’ 

q 


m/ \ 27r£ f dz 

Tit) = -= e / rit \ ^ ■ 

wi J Q{z)fx{z) 

-q 


(58) 


Note that P{t) and T{t) have the meaning of the local (i.e. defined at a point (0,t)) periods of the slowly modulated 
finite-band solution - see the definition (20) of the local wavenumbers and local frequencies. Thus, the solution at 
x = 0 considered for any fixed time t within the genus two region has a single local spatial period and thus, can be 
approximated by an appropriate periodic NLS solution in some vicinity of x = 0. The natural candidate for such 
an approximation is the Akhmediev breather (AB) (11), which is a limiting wave form of the two-gap NLS solution 
and has a single spatial period 2'ne/p (we note that AB (11) is also a single-parameter solution so its period also 
defines the amplitude). To this end we use the dependence P{t) as the period for the AB and compare the profiles 
of intensity in the quasiperiodic breather lattice observed in the numerical solution (see Figs. 5, 8 ) with the spatial 
profile of the AB (11) with the period 27re/p = P(t) and appropriately chosen phase. Such a comparison for t = 0.676 
(P = 0.154) is shown in Fig. 10b. One can see excellent agreement for the amplitudes, positions and detailed profiles 
of the main peaks (obviously some fitting of the AB phase was necessary). Remarkably, the agreement within the 
genus two region remains very good (and even improves with respect to the lower maxima) away from x = 0. Thus, 
the amplitude profile of the breather lattice in the bulk of the genus two region can be approximated by the modulus 
of the modulated “time-periodic AB” with the slowly varying spatial and temporal periods given by (58). 

As one can see, the spatial period P increases with time which is indicative of the tendency of the approximating 
AB breather towards the Peregrine soliton. Indeed, from the trajectory of the branch point a in the complex plane 
shown in Fig. 9b it is seen that near the upper boundary of the genus two region, (t « 1.25 - see Fig. 5), the imaginary 
part b of the moving branch point approaches the value q = 1 whereas a approaches zero. Since a = a 2 = —ai one 
can conclude that the spectral portrait of the solution approaches that of the Peregrine soliton (two double points 
placed at the endpoints -Piq of the basic branchcut). Indeed, the comparison of the wave form of the large oscillation 
observed in the numerical simulations at f = 1.16, x = 0 with the plot of the amplitude |'0| for the Peregrine soliton 
(14) shows excellent agreement - see Fig. 11b (we have also checked the agreement for the phase but do not show it 
on the plot). 

Finally, we study the behaviour of the maximum wave height as a function of time within the genus two region. 
The maximum wave hight in a (non-modulated) finite-band NLS solution with 0 < (/ < 2 can be found from the 
simple formula \'4’\m = where bj = Imja^] [79]. This formula is obvious for g = 0,1 (see ( 8 )). It can also 

be readily obtained for g = 2 in the particular case when bi = 62 (see e.g. [5]). To this end we plot the function 
1 -I- 2b{t) using the solution b{t) of the modulation equation (56) and compare the result with the values of IV'lm 
extracted from the numerical solution of the box problem for the NLS equation (1) with very small s = 1/60 (see 
Fig. 5b for the corresponding amplitude density plot) in the strip —0.2 < x < 0.2 containing at least three peaks 
in an e-neighbourhood of most points (0,t) within the g = 2 region. The comparison is presented in Fig. 12. One 
can see that the curve 1 -I- 2b{t) provides a very good approximation for the dependence of \ip\rait), which exhibits 
rapid growth within the two-phase interaction region, further supporting the proposed mechanism of the rogue wave 
formation. 
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FIG. 12: (Color online) Maximum amplitude \ip\m as function of time in the genus two interaction region. Solid line: the 
estimate 2b{t) + 1 for |V'|m constructed from the modulation solution (56) at a; = 0; Dots: the values of extracted from 
the numerical solution of the NLS box problem with q — I, L = 25/33, e = 1/60 (see Fig. 5b) within the strip —0.2 < x < 0.2. 


3. Long-time behaviour 

Of particular interest is the long-time behaviour of the solution to the semi-classical NLS box problem. Here we only 
present a hypothesis about the asymptotic structure of the solution based on the results of the numerical simulations, 
leaving a detailed analytical study to future work. 

The small-dispersion NLS evolution in the box problem (1), (38) with zero initial ‘velocity’ u has two key macroscopic 
features: (i) the oscillations for all times are confined to the spatial domain —L < x < L\ and (ii) the solution genus 
(the number of nonlinear oscillatory modes) grows with time. Both features are illustrated in the diagram in Fig. 6a. 
Our numerical simulations suggest that the pattern of the x-t plane splitting into the regions of different genera shown 
in Fig. 6a persists as time increases (we ran the computations up to about t = 6 for the box problem with e = 1/33, 
q = \, L = 25/33). Motivated by this observation, we put forward a plausible hypothesis that for t ^ 1 the solution 
genus g t, as long as t <C e~^. Then a pertinent question arises: what is the long-time asymptotic distribution of the 
spectral bands in the complex plane? In the KdV theory the consideration of the thermodynamic type infinite-genus 
limit of finite-band potentials [31] and of the associated Whitham equations [32] has lead to the kinetic description 
of a soliton gas [33], [34]. The focusing NLS counterpart of this theory would include the ‘breather gas’ description 
which is yet to be developed. The long-time asymptotics of the NLS box problem, thus, could provide insight into 
the properties of strong integrable NLS turbulence, which has recently become the subject of an active research (see 
[3] and [20] for the recent numerical and experimental results respectively). 



t 



X 


FIG. 13: (Golor online) Formation of higher order rogue waves in the long-time evolution in the NLS box problem with q — 1, 
L = 25/33, e = 1/33. a) Dependence of the wave amplitude \tp\ on t at a: = 0. One can see the presence of the rogue wave with 
1^1^ « 4.5 at t = 4.118; b) Profile of the higher-order rogue wave: numerical solution of the box problem at t = 4.100; 4.118 
(dashed and dotted-dashed lines respectively) near a; = 0 and second-order rational breather solution (59) (solid line). 
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In the higher-genus {g > 2) regions generated in the course of the evolution in the NLS box problem, there arises 
a possibility of the formation of higher-order rogue waves with the maximum height signihcantly exceeding that of 
the ‘regular’ rogue waves. Probably the simplest example of such a “super rogue wave” is the second-order rational 
breather solution of the NLS equation, which has the form [75], [15] 


'0 = 9 


1-4 


G + iH 
D 




where G, H and D are given by: 


3 9X^ 

64 liT ~T 



x^ 

~Y 


16 2 


f ^X^ +X'^ + -h 6V^T^ + , 

- 3X2 2X2 ^4 4X2r2 -H 2T^^ , 

X 2 y 2 9 J.4 ^2j,4 ^ ^ 

4 3 


(59) 


(60) 


Here X = cc/e, T = t/e. The maximum hight of the breather (60) is 5g. 

Our numerical simulations show that the higher order rogue wave indeed appear in the NLS box problem. In 
Fig. 13a the plot of the numerical solution for the amplitude ['0(0, f)| is presented for the box problem with e = 1/33, 
q = 1 and L = 25/33 in the interval 0 < f < 6 . One can see a very high peak with j'i/jm ~ 4.5 at about t = 4.1. 
The comparison of the wave profile at t = 4.118 with the second-order rational breather profile (60) is shown in 
Fig. 13b and demonstrates good agreement. This ‘super rogue wave’ can be interpreted as a result of the collision of 
two lower-order breathers shown in dashed line corresponding to the solution at t = 4.100, just prior the formation 
of the higher order rogue wave. A more accurate interpretation of this effect is the interaction of nonlinear modes 
within the modulated g-phase ( 5 -band) solution. In this regard it is worth noting that, while the local profile of 
this solution around a; = 0 is reasonably close to the rational breather, it is not an exact solitary wave (see also the 
comparisons with Akhmediev and Peregrine breathers in the 5 = 2 region). Importantly, the emergence of a single 
large amplitude oscillation within a multiphase solution (4) at a particular cc, t-point requires that a certain precise 
relationship between the phases s~^r]j is satisfied, and so is sensitive to changes of e. This sensitivity is expected to 
increase with growth of 5 . Thus, the exact prediction of the emergence of higher-order rogue waves in the regions with 
sufficiently large 5 is impractical and should be replaced with a statistical description within the general integrable 
turbulence theory even though the original formulation of the semi-classical NLS box problem is purely deterministic. 
This is in line with the proposition made in the beginning of this section that the long-time asymptotic behaviour 
of the semiclassical NLS should be generally described in statistical terms. In this connection we note that the 
statistical description of the long-time asymptotic solution of the small-dispersion KdV equation with deterministic 
initial conditions defined on the entire cc-axis was considered in [80], [81]. 


V. EFFECTS OF PERTURBATIONS ON THE SEMI-CLASSICAL EVOLUTION 

In Sections II - IV we have described an analytically tractable scenario of the rogue wave generation in the framework 
of the semi-classical focusing NLS equation with the inital data in the form of a real-valued rectangular potential (the 
‘box’). In practice (physical or numerical experiment) this idealised scenario may be affected by at least two factors: 
(i) the presence of a noise (physical or numerical); and (ii) higher-order physical effects (e.g. Raman scattering, 
saturable nonlinearity etc). The first factor is inevitably present in any physical system and, due to modulational 
instability, will impose natural restrictions on the admissible values of q and L characterising the initial box potential. 
The second factor generally destroys integrability of the NLS equation and thus, can affect the very existence of the 
multi-phase solutions and the corresponding semi-classical limits. While, obviously, the quantitative effect of both 
factors depends on their magnitudes, it is important to understand what qualitative changes may occur in the system 
due to their presence. While the detailed study of this important issue is beyond the scope of this paper, we present 
below some estimates and numerical simulations illustrating the effects of the noise and non-integrable perturbations 
on the solution of the NLS box problem. 

The effect of the external noise on the dispersive dam break flow evolution was briefly discussed at the end of 
Section III (see formula (37) for the estimate of the dispersive dam break flow lifetime due to the development of 
modulational instability of the external condensate (plane wave)). In the box problem involving the generation of 
two dispersive dam break flows, the presence of the noise will impose some restrictions on the admissible initial box 
parameters for which the semi-classical NLS description is valid in “practical terms”. E.g. if the box is too wide or 
too tall, the noise perturbations of the condensate in the central part of the box will have enough time to develop 
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before the harmonic edges of the counter-propagating dispersive dam break flows will meet at a; = 0 and saturate 
the instability. In Fig. 14 the amplitude density plot is presented for the NLS evolution of the initial box with our 
standard parameters q — 1^ L = 25/33 but with e = 0.01 in the NLS equation (1), which is significantly smaller than 
our usual value e = 1/33. One can see an extra oscillatory structure forming in the central part of the box, inside the 
region g = 0. Apparently, a related phenomenon was observed in the numerical simulations in [50] where it was aptly 
named “the beard”. In [50] the beard phenomenon was ascribed to non-analyticity of the initial data. Our numerical 
experiments with the box initial conditions suggest that the origin of the beard phenomenon lies in the development 
of the modulational instability due to the presence of a numerical noise. 




FIG. 14: (Color online) “Growing the beard”: the effect of the noise on the semi-dassical box evolution. Numerical simulation 
of the NLS box problem with g = 1, L = 25/33, e = 0.01. a) density plot for b) amplitude profiles for different t. 


We use the estimate (37) for the lifetime tm of the dispersive dam break flow to obtain an estimate for the “critical” 
parameters L^, Qc of the initial box with respect to the beard formation phenomenon. For that, we use the balance 
relation = toi where to = L/2\/2q is the time at which the collision of two dispersive dam break flow occurs at 
a: = 0 (see (41)) so that the modulational instability of the plane wave in the genus zero region is saturated by the 
formation of breather lattices. As a result we obtain the critical value of the initial “mass” A = qL'. 

Ac = (qL)c = V2eln ^ , (61) 

0 

where S is the typical amplitude of the noise. The boxes with A > A,, will “grow” the beard. 

The higher order physical effects are described by additional/modified terms in the NLS equation. In most cases 
these modifications lead to the loss of integrability and hence, the 1ST and the semi-classical analysis via the Riemann- 
Hilbert steepest descent approach are no longer available. Still, periodic solutions and the Whitham equations can be 
derived, so the formal analytical description of dispersive dam breaks flows is possible. It is interesting to see whether 
the interaction pattern for counter-propagating dam break flows will persist despite non-integrability of the problem. 



FIG. 15: (Color online) Effect of nonlinear saturation on the formation of a breather lattice in the small-dispersion box problem 
(cf. Figs. 5, 8 ). Numerical solution of the box problem {q = 1, L = 25/33) for the saturable NLS equation (62) with 7 = 0.1, 
e = 1/30. a) Density plot for \ij}{x,t)\', b) Spatial profiles of the amplitude ['(/’(s^)! a-t different moments of time. 
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To be specific, we consider the NLS equation with saturable nonlinearity, which arises in the description of light 
propagation in some media with highly nonlinear optical properties (see e.g. [82]) 

^ ° 

where 7 is the coefficient characterising the strength of the saturation effect. For 7^1 equation (62) is a perturbed 
NLS equation ( 1 ). 

In Fig. 15 the results of the numerical simulations of the box problem with q = 1, L = 25/33 for the saturable NLS 
(62) with 7 = 0.1, e = 1/30 are presented. We have chosen a relatively large value of the parameter 7 to elucidate the 
qualitative effects of saturable nonlinearity. One can draw some immediate conclusions from the simulations shown 
in Fig. 15: (i) the initial evolution for the saturable case is qualitatively similar to that in the pure, cubic NLS case: 
one can clearly see the formation of two dispersive dam break flows (cf. Fig. 5); (ii) quite remarkably, the qualitative 
agreement with the cubic case stretches beyond the evolution of the modulated periodic solutions. Indeed, one can 
see that the interaction of two dispersive dam break flows leads to the formation of the two-phase region dominated 
by the breather lattices despite non-integrability of the saturable NLS (62) (we note that a similar persistence of the 
two-phase pattern despite non-integrability of the governing equation was observed in the numerical simulations of 
the DSW interaction for the defocusing saturable NLS [83] and in the analysis of DSWs in viscous fluid conduits [84]). 
One should also note that the amplitudes of the breathers generated in the saturable NLS are noticeably smaller than 
in the cubic nonlinearity case (cf. Fig. 15b and Fig. 8 ); (hi) the evolution beyond the two-phase region is qualitatively 
and quantitatively very different compared to the integrable case. 


VI. CONCLUSIONS 

In this work, we have proposed a novel mechanism of the rogue wave formation described in the framework of the 
focusing nonlinear Schrodinger (NLS) equation with small dispersion parameter. The key role in our construction is 
played by the dispersive focusing dam break flow — a DSW-like nonlinear wave train regularising an initial sharp 
transition between the uniform plane wave and the zero-intensity, vacuum state. We have considered the NLS 
evolution of a square profile (a “box”) giving rise to two such counter-propagating dispersive dam break flows, whose 
interaction has been shown to result in the emergence of a modulated two-phase large-amplitude breather lattice 
closely approximated by a sequence of Akhmediev and Peregrine breathers within certain time-space domain. We 
have used a combination of the nonlinear modulation (Whitham) theory and elements of the steepest descent method 
for the Riemann-Hilbert problem associated with the semi-classical NLS equation, to construct the exact modulation 
solution describing the two-phase interaction in the box problem, and predict the parameters of the emerging rogue 
waves. Our semi-classical analytical results are shown to be in excellent agreement with direct numerical simulations of 
the small-dispersion NLS box problem. We also show that the proposed rogue wave generation mechanism is different, 
both physically and mathematically, from the generation of the Peregrine breathers during the regularisation of the 
generic gradient catastrophe in the semi-classical NLS equation with analytic, bell-shaped initial data considered in 
[12, 44-46]. 

Our numerical simulations of the NLS box problem for longer times suggest that the evolution beyond the two- 
phase interaction dynamics leads to the generation of the regions filled with quasi-periodic waves with the number 
of interacting nonlinear modes (phases) increasing with time at each point within the spatial location of the initial 
box potential. We put forward a hypothesis that for 1 ^ t <C £~^ the number of oscillating phases (the genus of 
the solution) grows as g ^ t. It is argued that such a complex multi-phase wave structure would require a statistical 
description similar to that constructed in [31-33] for the KdV soliton gas/soliton turbulence. Such a description is 
yet to be developed and could provide an important insight into properties of the NLS integrable turbulence. 

Finally, we numerically considered effects of small perturbations on the qualitative structure of the small-dispersion 
NLS box problem solution. We first looked at the effect of noise, inevitably present in any physically (and numerically) 
realistic setting. This effect becomes essential for small values of the dispersive parameter e and is manifested in the 
generation of an additional oscillatory structure - the ‘beard’ - in the central part of the box. This structure is not 
captured by the semi-classical NLS box problem solution. Secondly, we performed numerical simulation of the box 
problem for the small-dispersion NLS equation with saturable nonlinearity, a perturbed version of the cubic NLS 
equation. Although the weakly saturated nonlinearity introduces a non-integrable perturbation, our simulations show 
that, remarkably, this does not destroy the qualitative structure of the breather lattice even for relatively large values 
of the saturation parameter, although the amplitude of the breathers is noticeably smaller than in the unperturbed, 
cubic, case. 

The proposed mechanism of the rogue wave formation can be realised in fibre optics experiments. The obtained 
analytical solutions for the interaction of dispersive dam break flows can also find applications in oceanography and 
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BEC dynamics (see the physical estimates in [12] and [44] showing the relevance of the small-dispersion focusing NLS 
to these two areas). 

The general conclusion to be drawn from this work is that the semi-classical NLS equation provides a powerful 
analytical framework for the description of physically important effects related to the rogue wave formation and 
transition to integrable turbulence. 
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Appendix A: Elements of the Riemann-Hilbert problem analysis for the semiclassical focusing NLS equation 

with the box-like initial data. 


1 . 0 -function 


The ID NLS equation with cubic nonlinearity considered in this paper is an integrable equation [16], which can 
be solved through the inverse scattering method. It was observed in [76] that the inverse scattering transform (1ST), 
which is in the core of the method, can be written as a (multiplicative) matrix Riemann-Hilbert Problem (RHP). 
In the study of the semiclassical (zero dispersion) limit of the NLS we are naturally interested in the corresponding 
e —0 limit of the 1ST. The recently developed nonlinear steepest descent method of Deift and Zhou [43] is a very 
effective tool for small or large parameter asymptotics of matrix RHPs. The nonlinear steepest descent asymptotics 
of the semiclassical NLS eq. (1) was first obtained in [71] (pure soliton case) and [64] (solitons and radiation). 
The central object of this analysis is the so-called g-function g(A) = g(A;a;,t), which is an analytic (and Schwarz- 
symmetrical) function on the Riemann surface Tg{x,t), whose branch points aj = aj{x,t) are determined by the 
Whitham equations. The framework of this paper does not provide enough room to properly describe g, we only 
mention that g(A;a;,t) satisfies the following jump conditions 

fl+(A)+g-(A) = 6»fc(A)-f 7?fc, fc = 0,1,...,5, (Al) 

where g± are the values of the g-function at the opposite sides of the oriented branchcut jk between ak and at', rik 
are some real constants (in A), and rjo = 0. Here 

6»fc (A; x, <) = /fc (A) -f 2tA^ -h 2xA, (A2) 

where the functions fkiz) contain the information about scattering data for a particular initial condition (potential). 
For the box potential [22] /o, /i and /2 are given by (52). 

In general, g may have additional constant jumps along some gaps connecting the neighboring bands, but in the case 
of the box potential there are no such jumps. Thus, the function g(A) is analytic everywhere except the branchcuts 
7 /c. In particular, it is analytic at A = oo. 

By the well known Sokhotski-Plemelj formula. 


2g(A) 


m) yi VidC f 0{OdC 


(A3) 


where the (Schwarz-symmetrical) radical Tl{X) = Ti-g{X) is given by (5), 7 j denotes the clockwise oriented loop around 
7 j and 6{X) is defined as 0k{X) on %. We assume the loops jj do not intersect each other and A is outside of any 
loop. Some mathematical details of the RHP derivations can be found, for example, in [64], [78]. Here we proceed 
with the brief presentation of the results relevant to the analysis of the semi-classical NLS box problem. 

The requirement that g(z) is analytic at 2 = oo leads to the following g real equations on gj: 
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27rf 



C^g(C) 

^(C) 


dC 



no 


dC = 0 


fc = 0,l,--- ,5-1- 


(A4) 
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To shorten the notation, we will proceed by considering the case of genus g = 2, with the understanding that all 
the following calculations can be readily generalised to an arbitrary genus 5 S N. Simple linear algebra shows that 
conditions {A4) combined with (A3) yield 


where 


200 ) = 




£ 

dc 

Jjl 

no 

£ 

dC 

*72 

no 

£ 

eiOdc 

Jujj 

no 


£ ^ 

^72 no 
r ceiQdc 

fujj no 


dC 


J 71 

ic-ono 

r e(Odc 

fuj, ic-ono 


D = 


dC 

(£ 

CdC 

no 

J71 

no 

dC 

<£ 

CdC 

no 

J72 

no 


The details of the calculations can be found in 


(A5) 


(A 6 ) 


2. RHP and modulation solution 


The Whitham modulation theory is concerned with the spatiotemporal evolution of the branch points aj 
through which the physical parameters of the solution (amplitude, wavenumber, frequency etc) are expressed. The 
functions aj(x,t), ckjix^t) are solutions of the Whitham modulation equations (17) with appropriate initial or bound¬ 
ary conditions derived from a given NTS initial-value problem. The RHP approach yields these dependencies directly, 
by-passing the procedure of the integration of the Whitham modulation equations. We stress that the modulation 
solution aj{x,t), aj(x,t) is only part of the full RHP solution, and finding the full solution could be a rather de¬ 
manding mathematical task. Thus, as long as the evolution of the branchpoints is concerned, one can take advantage 
of the appropriate part of the full RHP analysis for the derivation of the dependencies aj{x,t) and then verify the 
formal result by checking its consistency with the Whitham modulation equations and the corresponding initial or 
boundary conditions. A detailed analysis of interrelations between the Whitham modulation theory for the focusing 
NTS equation and the RHP approach can be found in [65]. 

It is shown in [77] that the equations for the moving (non-constant) branchpoints aj{x,t) follow from the condition 

K{aj)=0, where j = 1,2 and c.c. (A7) 

(In the case of box potential Oq = iq is a, fixed branch point). We note that equation (A7) for each particular aj is 
equivalent to [77] 

, 0 . (A8) 


Thus, on the solutions of the modulation equations (A7), 


d ^ d d ^ d 
dx^ dt^ dt^ 

In view of (A2), (A6) and (A9), the modulation solution (A7) can be written in the form 


(A9) 


X + 


Ktjaj) ^ ^ Ko{aj) 
Kx{oij) Kx{oij) 


and 


J = l,2, 


(AlO) 


where Kx,Kt denote partial derivatives of K, and Kq{z) is obtained from K{z) by replacing x,t with zero in (A2), 
(A6). 

Comparison of (AlO) and the solution of the Whitham equations in the generalised hodograph form (42) suggests 
the identification 


Vj{a,a) 


Ktjaj) 

Kx{aj) ’ 


Wj{a, a) 


Ko{aj) 
Kx {o-j) 


i = i,2. 


(All) 


Direct substitution shows that Tsarev’s equations (43) are indeed satisfied by (All). Thus, the branchpoints aj 
defined by (A7) satisfy the Whitham equations (17). The NTS initial conditions (the box potential (38) in our case) 
enter the modulation via the functions /fc(A) (52) defining the jump functions Ok (A2). The modulation solution (A7) 
(or, equivalently, (AlO)) constructed in such way automatically satisfies all the necessary matching conditions for a^, 
aj at the second breaking curve. We shall later verify the matching explicitly for x = 0. 
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3. Phases and characteristic speeds 


Our aim in this subsection is to obtain explicit expression (18) and (49) for the characteristic speeds Vj(a,a) and 
phases Tj(a,a) respectively. For that, we introduce the basis holomorphic differentials on the Riemann surface F by 


drij = 


PjW 

TZ{a, a, A) 


dX, 


J = l,2, 


(A12) 


where Pj{z) = >fj,iA + and the coefficients >Cj^k are found from the normalisation (7). Then (7) and (A6) imply 


f >(l,2 >(2,2\ _ jj-l 

x:2,ij ~ 


1^1 \ 4 7^ 4 7^ / 


(A13) 


Note that the coefficients Xij determine the wavenumbers and frequencies of the multi-phase solution, see (6). Making 
appropriate linear combinations of the first two columns in the determinant representation (A6) for K(X), we obtain 


K{X) 




1 
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piiOeiQdc 
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0 
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4 ic-ono 

HQdC 


ic-ono 


(A14) 


which, together with (A3) and (A5) implies that 


Vi = 


U7j 


mpAQdC 

no 


= Ani [Pjt + ~ > J = 1) 2, 


where (cf. (6)) 


Po - 




(A15) 


(A16) 


and Tj is given by (49). Here the residue theory was used to calculate (A15), (A16), (49). Now, using (A16), (A14) 
it is not difficult to show that the determinant form (AlO) of the modulation solution is equivalent to the hodograph 
solution (46) obtained in terms of the phase T^. We recall that for the box potential the functions fk{X) in (A16) are 
given by (52). 

Equations (A14), (A9) and (A2) yield convenient formulae 


K^{X;x,t) = i 


dC. 


Kt{X;x,t) = -2\D\J2[nn^ 

V 


N 


J2ak 


Lfc=l 


i=i 


+ 


(C - A)77(C) 
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dC, 

dC, 


(A17) 


so that the characteristic speeds (All) are 


) 


U-i 




+ 


(A18) 


Efc=i ^k,i (c-4mo^^ 

As it was mentioned in Section II, this expression can be readily generalized to the case of arbitrary genus g, see (21). 


4. Modulation solution for the box problem (g = 2, x = 0) 

We now look closer at the modulation solution (AlO) for the box potential. Specifically, we consider behaviour of 
the solution at a; = 0 in the genus two region. First, we recall that, for the box potential (38) we have ag = iq so 
7^(A) = i?(A)z^(A), where i?(A), iy{X) are given by (51). 

In the case x = 0 we have symmetry 02 = —di. Then i?(A) = 0 (A^ — a^)(A^ — d^), where a. = 02 - 
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By taking linear combinations of the first two rows in (A6), we can make 72 ± 71 to be the countours of integration 
for the first and for the second row of (A6). Then 


K{\)=4. 


dC 
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r dC 
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27r2 J72 (C-A)fl(C) 


(A19) 


where the contours 72 ± 71 were deformed to [—iq,iq] and iM. \ [—iq,iq] respectively, and the subscript ‘+' in 
indicates the limiting value of on the positive side of the branchcut {—iq, iq\. Note that the function i?(A) is even 
on zK whereas z^+(A) is even on [—iq,iq\ but odd on zR \ [—iq,iq\. Thus, 
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so the complex modulation equation becomes 
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(A20) 


We denote a = a + z6 and make the change of variable ^ = zz to represent (A20) in the form (56) containing only 
integrals over intervals of real line. 

As a by-product of our calculation we derive explicit expressions at x = 0 for the coefficients of the holomorphic 
differential (A12). These coefficients determine, via (6), the local wavenumbers and frequencies of the modulated 
multi-phase wave. 

For g = 2 the values >Ci j are defined by the general formula (A13) following from the normalisation conditions (7). 
The crucial observation is that the determinant |iA| in (A13) is the (3,3) minor of the main determinant K{z) (see 
(A6)). Then, at x = 0 we take advantage of the representation (A20) for K{a) to obtain: 
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9] ^(C)^+(C) 
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CdC 

-iqyq] RiCMC) ' 


Now, using (A13) and taking into account the symmetry 02 = —ai we obtain 


X = 0 : >!:i.i(a, a) 


2/ 
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Jl-zqyq] RiCMO 



(A21) 


(A22) 


where, we recall, a = a^- Introducing the change of variable C, = iz we represent (A22) in the form (54) containing 
only integrals along the real axis, 


Appendix B: Numerical method 


Here we only present a brief description of the numerical method used in this paper for solving the small dispersion 
NTS equation (1), leaving details to a separate publication. We first scale the time variable in (1) through t = 2t and 
write = u + iv, where u and v are real-valued functions. Then, Eq.(l) can be written as the following system of 
equations: 


Ur = -eVxx - [u +v ) v , 

e ^ ' 

Vr = eUxx + - ( u ^ + u . 

e ^ ' 


(Bl) 
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The time derivatives Ut and Vr in (Bl) are found by the 4th-order Adams-Bashforth-Moulton (ABM) predictor- 
corrector method [85]. The first four time steps are solved by another method (e.g. 4th-order Runge-Kutta) since 
the ABM method needs four initial values to be started. The spatial derivatives v^x and Uxx are calculated using 
a pseudo-spectral derivative approximation without any filtering. Thus, the resulting algorithm is simple, totally 
explicit and can provide long-time numerical evaluation without generating numerical artifacts for reasonably small 
values of e, e.g. e = 1/33. The stability region for the ABM method is narrower than that for the traditional Fourier 
split-step method, widely used for solving the defocusing NTS equation. However, this latter method can easily yield 
wrong results when using small values of e. 

The above algorithm has been tested with both discontinous and smoothed out initial data. We have verified that 
the smoothing clearly preserves the structure of the solution with the advantage of a better control, as the value of e 
is decreased, over the round-off errors than in the discontinuous data case (the derivatives in the regions containing 
discontinuity are difficult to approximate numerically). As a result, one can avoid the (numerically induced) effects of 
modulational instability which are unavoidable with the discontinuous initial data. The only drawback of smoothing 
the initial data is that the formation of the edge soliton is delayed and, as a result, it has a different phase compared 
to the discontinuous case. 
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